%% Likelihood Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function nval = lik_fcn(prmtr,T,y,mstarU,mstarL,vol_bdate,growth_bdate)

           prmtr = trans(prmtr);

% Number of states to keep track of            
           no_st=8;               % max number of periods to be considered; need to adjust lik_fnc to allow more than 8 
           dmnsion=3^no_st;       % max number of cases to be considered 
    
% Parameters, where regimes: 0 normal; 1 L-shaped recession; 2 U-shaped recession 

           p01 = prmtr(1);      %Pr[St=1|St-1=0]
           p02 = prmtr(2);      %Pr[St=2|St-1=0]
           p11 = prmtr(3);      %Pr[St=1|St-1=1]
           p22 = prmtr(4);      %Pr[St=2|St-1=2]      
           mu0 = prmtr(5);      % expansion growth 
           mu1 = prmtr(6);      % L-shaped recession effect 
           mu2 = prmtr(7);      % U-shaped recession effect
           omega1 = 0;           % no bounceback effect for L-shaped recession         
           omega2 = -mu2/mstarU; % bounceback effect for U-shaped recession 	   
           delta = prmtr(8);     % change in trend growth 
           sigma2_pre = prmtr(9)^2; % volatility before break
           sigma2_post = prmtr(10)^2;  % volatility after break
		   
% A MATRIX OF TRANSITION PROBABILITIES 

           PI = [1-p01 - p02, 1-p11, 1-p22;
                    p01, p11, 0;
                    p02, 0, p22];
              
% INITIALIZING THE FILTER WITH STEADY-STATE PROBABILITIES 

           A = [eye(3)-PI; ones(1,3)];
           en=[zeros(3,1);1];
           prob__t = inv(A'*A)*A'*en;  %PR[S_t=0]|PR[S_t=1]|PR[S_t=2],3x1 steady-state probs
                             
                                                                                              
           pr_trf0 = reshape (PI, [numel(PI),1]);
           pr_trf1=[pr_trf0;pr_trf0;pr_trf0];
           pr_trf2=[pr_trf1;pr_trf1;pr_trf1];
           pr_trf3=[pr_trf2;pr_trf2;pr_trf2];
           pr_trf4=[pr_trf3;pr_trf3;pr_trf3];
           pr_trf5=[pr_trf4;pr_trf4;pr_trf4];
           pr_trf=[pr_trf5;pr_trf5;pr_trf5];
                      

           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf0;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf1;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf2;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf3;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf4;
           prob__t = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]) .* pr_trf5;
           prob__  = reshape ([prob__t, prob__t, prob__t]', [numel([prob__t, prob__t, prob__t]'),1]); % unconditional prob for 3^8 cases
    
            sU = [0;0;1];

            s0 = kron(ones(3^7,1),(kron(sU,ones(3^0,1))));
            s1 = kron(ones(3^6,1),(kron(sU,ones(3^1,1))));
            s2 = kron(ones(3^5,1),(kron(sU,ones(3^2,1))));
            s3 = kron(ones(3^4,1),(kron(sU,ones(3^3,1))));
            s4 = kron(ones(3^3,1),(kron(sU,ones(3^4,1))));
            s5 = kron(ones(3^2,1),(kron(sU,ones(3^5,1))));
            s6 = kron(ones(3^1,1),(kron(sU,ones(3^6,1))));
            s7 = kron(ones(3^0,1),(kron(sU,ones(3^7,1))));

            sUmat = [s7,s6,s5,s4,s3,s2,s1,s0];

            sL = [0;1;0];

            s0 = kron(ones(3^7,1),(kron(sL,ones(3^0,1))));
            s1 = kron(ones(3^6,1),(kron(sL,ones(3^1,1))));
            s2 = kron(ones(3^5,1),(kron(sL,ones(3^2,1))));
            s3 = kron(ones(3^4,1),(kron(sL,ones(3^3,1))));
            s4 = kron(ones(3^3,1),(kron(sL,ones(3^4,1))));
            s5 = kron(ones(3^2,1),(kron(sL,ones(3^5,1))));
            s6 = kron(ones(3^1,1),(kron(sL,ones(3^6,1))));
            s7 = kron(ones(3^0,1),(kron(sL,ones(3^7,1))));

            sLmat = [s7,s6,s5,s4,s3,s2,s1,s0];
      
        
    val=0;
    
    d = 0; dd=0;

    j_iter = 1;
    while j_iter <= T % Loop over time periods
            
            if j_iter > growth_bdate        
                d = 1;    
            end  
            
            if j_iter > vol_bdate        
                dd = 1;    
            end  
            
         
        f_cast1 = y(j_iter)*ones(dmnsion,1) - ones(dmnsion,1)*mu0 - d*ones(dmnsion,1)*delta - sLmat(:,8)*mu1 - omega1*sum(sLmat(:,8-mstarL:7)')' - sUmat(:,8)*mu2 - omega2*sum(sUmat(:,8-mstarU:7)')';
                                             
        var_L = (1-dd)*sigma2_pre*ones(dmnsion,1) + dd*sigma2_post*ones(dmnsion,1);                                    

        prob_dd=pr_trf .* prob__;
                    % Pr[S_t,S_{t-1},.....,S_{t-7} | I(t-1)]

        pr_vl=(1./sqrt(2.*pi.*var_L)).*exp(-0.5*f_cast1.*f_cast1./var_L).*prob_dd;
                                                                 % 3^8 x 1 Joint density of y_t,S_t,..,S_{t-7} given past information

        pr_val=sum(pr_vl');         % f(y_t|I(t-1)), density of y_t given past information:  This is weighted average of 3^8 conditional densities

        lik=log(pr_val);

        pro_=pr_vl/pr_val; % Pr[S_t,S_{t-1},...,S_{t-7} | I(t-1),y_t]
                      % Updating of prob. with new information, y_t   


        prob__t=pro_(1:dmnsion/3,1)+pro_(dmnsion/3+1:2*dmnsion/3,1)+pro_(2*dmnsion/3+1:dmnsion,1);
                       % Integrate out S_{t-7}: then you get Pr[S_t, S_{t-1},....,S_{t-6}| Y_t]  
                       %  3^7 x 1

        prob__=reshape([prob__t,prob__t,prob__t]', [numel([prob__t,prob__t,prob__t]'),1]);  
            
        val = val+lik;

   j_iter = j_iter + 1;
   end % End of loop over time periods

   nval = -val; % Negative log likelihood value
    
end % End of likelihood function